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CN Abstract 

xyy In this paper, we describe a multiscale model based on magneto-static traps of 

Q neutral atoms or ion traps. The idea is to levitate a magnetic spinning top in the 

air repelled by a base magnet. 

For such a problem, we have to deal with different time and spatial scales and we 
propose a novel splitting method for solving the levitron problem, see [1]. 

We focus on the multiscale problem, which we obtain by coupling the kinetic T 
7— I and the potential U part of our equation. The kinetic and potential parts, can be 

seen as generators of flows, see [2]. 

The main problem is based on the accurate computation of the Hamiltonian 



equation and we propose a novel higher order splitting scheme to obtain stable states 
near the relative equilibrium. To improve the splitting scheme we apply a novel 



q method so called MPE (multiproduct expansion method), see [3], which include 

CN higher order extrapolation schemes. 

In numerical studies, we discuss the stability near this relative equilibrium with 
.£h our improved time-integrators. Best results are obtained by iterative and extrapo- 

lated Verlet schemes in comparison to higher order explicit Runge-Kutta schemes. 
Experiments are applied to a magnetic top in an axisymmetric magnetic field (i.e. 
the Levitron) and we discuss the future applications to quantum computations. 
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1 Introduction 



We are motivated to simulate a Levitron, which is a magnetic spinning top and 
can levitate in a magnetic field. The main problem of such a nonlinear problem 
is to achieve a stability for the calculation of the critical splint rate. While the 
stability of Levitrons are discussed in the work of [4] and their dynamics in 
[5] , we concentrate on improving the standard time-integrator schemes for the 
Hamiltonian systems. It is important to derive stable numerical schemes with 
high accuracy to compute the non-dissipative equation of motions, which are 
at least higher order symplectic integrators. Here, we apply geometric integra- 
tors based on the Stromer- Verlet method with extrapolation methods, see [3] . 
While we have symplectic schemes, we preserve the underlying physics of our 
Levitron, i.e. reversibility, symplecticity, volume preservation and conservation 
of the first integrals, see [6]. 

For the numerical studies, we propose novel splitting schemes and analyze 
their behavior. We deal with a standard Verlet integrator and improve its 
accuracy with iterative and extrapolation ideas. Such a Hamiltonian splitting 
method, can be seen as geometric integrator and saves computational time 
while decoupling the full equation system. 

The paper is organized as follows. A mathematical model based on a multiscale 
problem of the Levitron is introduced in Section 2. The splitting method is 
used as a solver method to decouple the multiscale equations to more simpler 
equations given in the kinetics and the potential of the models is described 
in Section 3. The improvement of the splitting strategies based on extrap- 
olation schemes are discussed in Section 4. The numerical experiments and 
their description of our used methods are described in Section 5. Finally the 
conclusions and on overview for our next works are discussed in Section 6. 



2 Mathematical model 



The Levitron is described on the base of rigid body theory. With the con- 
vention of Goldstein [7] for the Euler angles the angular velocity is along 
the z-axis of the system, ug along the line of nodes and Coy, along the z'-axis. 
Finally the kinetic energy can be written as 



T = 



m(x 2 + f + z 2 ) + A{9 2 + 2 sin 2 9) + C(ip + <p cos 9f 



(1) 



The potential energy U is given by the sum of the gravitational energy and the 
interaction potential of the Levitron in the magnetic field of the base plate: 
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U = mgz — u,(smib sin# h cos ib sin 9 hcos# — ) (2) 

x y z 



with jj, as the magnetic moment of the top and <3> the magneto-static potential. 
Following Gans [5] we uses the potential of a ring dipole as approximation 
for a magnetized plane with a centered unmagnetized hole. Furthermore we 
introduced a nondimensionalization for the variables and the magneto-static 
potential: 



# = 



(1 + Z 2 ) 3 /2 



2 3 (2Z 2 -3)Z 

1 + M(l + Z 2 )7/2 



(3) 



Lengths were scaled by the radius R of the base plane, mass were measured in 



units of m and energy in units of mgh. Therefore the one time unit is JR/g. 



Knowing the kinetic and the potential energy, we can formulate the the La- 
grangian as: 



C=T(x 2 ,y 2 ,z 2 )-U(x,y,z) 

= ^ [m(x 2 + y 2 + z 2 ) + A(9 2 + 2 sin 2 9) + C(ip + <p cos 9) 2 
+ fi [sin 9 sin ipB x + sin 9 cos ipB y + cos 9B Z ] — mgz 



(4) 



Furthermore, the Hamiltonian H can be calculated as 



2m 



pi 



2A 2C 



P% . {P^-P-4> 



COS I 



2 A sin 2 9 

dtp dip d<p 

sin # sin-?/;— — V sin 9 cos ip— — V cost*—— 
ox ay dz 



+ mgz . 



(5) 



where 



q =(x;y;z;9;ip;<j)) 

_ (Px.Py^ Pz_,Pe.P± _ P<t> cos 9 -p^ cos 2 9 _ p<f, -p^ cosfl " 
" m'm'm'I'C Asin 2 # ~' Asm 2 9 



(6) 
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and p is given as: 





dC 


Px = 


dx 




dC 


Py = 


dy 




dC 


Pz = 






dC 


Pe = 


~m 




dC 


= 


dip 




dC 


P4> = 





mx 



= my 



mz 



A9 

C(ip + cos 9) 

A4> sin 2 9 + C(ip + cos 9) cos 9 



(7) 
(8) 

(9) 
(10) 

(11) 
(12) 



The Hamiltonian H is calculated as 



H =?p- C 



1 ( 2 . 2 . 2\ , Pe , Pi , (P<p-Pip 



2m 



COS I 



2,4 sin 2 

eta dcp dcp 

sin 6 1 sin ip — — h sin # cos -^t; — h cos#— — 
ax ay az 



(13) 



In the next section, we discuss the time-integrator methods to solve our dif- 
ferential equations. 



3 Splitting Methods 



The evolution of the dynamical variable u(q, p) (including q and p themselves) 
is given by the Poisson bracket, 



a / \ f9u dH du dH\ . 



A and B are Lie operators, or vector fields 



A = 



dH d 



B = 



dH d 



(14) 



dp <9q <9q dp 

The transfer to the operators are given in the following description. 



(15) 
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The exponential operators e AtA and e AtB are then just shift operators, with 
% (At) is a symmetric second order splitting method: 

r 2 yv (At) = e( A '/ 2 ) B e AtA e^ B . (16) 
and corresponds to the velocity form of the Verlet algorithm (VV). 
Further the splitting scheme: 

T 2 ,py(At) = e( At ^ A e AtB e( At ^ A . (17) 

and corresponds to the position- form of the Verlet algorithm (PV). 

In the literature, see [6], they are also known as symplectic splitting methods, 

see Figure 1. 

r C(t" +1/2 ) ^ 0(^1/2) 

B/2 B/2 



c(f) A/2 c( ,nj 

SE1 SE2 



c(t" +1 ) c(f +1 ) 
fB/2 



B 



B/2 I / 



c(t") A/2 c(f) 
Stoermer- Verlet (PV) Stoermer- Verlet (VV) 

Fig. 1. Symplectic Splitting Methods. 

The symplectic Stormer- Verlet or leap-frog algorithm in the notation 72, (At) 
SE2(At/2) o SEl(At/2) is given in the following algorithm 1. 

Algorithm 1 We start with (qo,Po)* = (q(t n ), p(t n ))*-' 



1 f)¥f f) 

(qi, Pi)* = e A ^(q , po)' = (/ - 2 At £ — (p., q*)^)(qo, Po)<, (18) 



(q 2 , v 2 )* = e^( qi , Vl )* = (/ + At £ — (p„ q,)^)(qi, Pi)', (19) 



i f)fJ f) 

(03, V 3 )* = e At / 2B (q 2 , P2 )*=(/- -At^:— ( Pt , qi )— )(q 2 , P2 )*. (20) 



ylnd t/ie substitution is given the algorithm for one time-step n — > n + 1 and 
we obtain the solution in t n+1 : 
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(q(r+ 1 ),v(r+ 1 ))* = (q 3! V3) t . 



For studying the delicate higher accurate stability of the Levitron, it is neces- 
sary to improve the standard Stormer-Verlet schemes, which are only second 
order scheme, but can be improve to a more accurate higher order scheme. 
In the following, we discuss the extrapolation idea with respect to the basic 
Stormer-Verlet algorithms. 



4 Improvement of the Splitting schemes with Extrapolation meth- 
ods 



In the following, we discuss the different splitting schemes, that are based on 
the Strang-splitting scheme, see [8] and extrapolated with the so called MPE 
method, see [3] for time-dependent problems. In the first part we present the 
linear version, while in the second part we embed an iterative scheme to derive 
a nonlinear version. 

The solution to the differential equation (14) can be formally written as 



/ rt+At \ 

u(t + At) = T(expy A(s)dsju(t), (21) 
where we assume general time-dependent operators A(t) — A(t) + B(t). 

Following by Suzuki [9] , we have a forward time derivative operator, also called 
super-operator: 

B = l < 22 > 

such that for any two time-dependent functions F(t) and G(t), 

F(t)e AtD G(t) = F(t + At)G(t). (23) 

If F(t) = 1, we have 

le AtD G(t) = e AtD G(t) = G(t). (24) 

By comparing with Trotters formula we can apply the Suzuki's decomposition 
of the time-ordered exponential and obtain: 

/ rt+At \ 

T( exp J A(s)ds) = exp[At{A{t) + £>)]. (25) 
Thus time-ordering can be achieve by splitting an additional operator D. 
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With such a scheme, we can transforms in any existing splitting algorithms 
into integrators of non-autonomous equations. 



Corrolar 2 We achieve the following important second order symmetric split- 
ting scheme for the Stormer-Verlet (PV) scheme, see also Figure 1: 

T 2 (At) = e^ AtD e AtA(t) e^ AtD = e 5 AtjB( ' + i Ai) e A * A (*+2 A *) e i A * B (*+3 A *), (26) 

which is the second-order scheme with the assumption of the commutation of 
the A and B between the D operator, see also [10]. 

Proof 3 For the second order algorithm, we apply the Strang- splitting scheme 
for the three operators A(t), B(t), D and we have assumed: [A(t),D] = 0, 
[B(t),D] = 0. 



72 (Ai) = e^ AtD e 5 A * B (* + l A * D ) e A * A (*) e l Ats W e^ AtD 

= e §AtB(t) e |A*D e AtA(t+±At) Q \AtD Q ±AtB(t+\At) Q \AtD 
_ e |AtB(t) Q \AtD e &tA{t+\At) e ±AtD e \AtB{t+\At) Q \AtD 
= e \MB{t+\ At) e AtA(t+\At) e \AtB(t+\At) ^ ^7) 

where we have applied the commutativity with the D operator and the shift 
with the forward time derivative operators. 

Remark 4 Every occurrence of the operator e diAtD , from right to left, updates 
the current time t to t + diAt. Ift is the time at the start of the algorithm, then 
after the first occurrence o/e2 AtD ; time is t + \At. After the second e2 A ' D ; 
time is t+At. Thus the leftmost ea AtD is not without effect, it correctly updates 
the time for the next iteration, see also [9]. 

Thus the iterations ofT^At) implicitly imply 



T 2 2 (At/2) = el^+f At ) e § At -4(*+3 A *), (28) 
by inserting equation (26), we obtained: 

7 2 2 (At/2) (29) 

= e \AtB{t+\At) e |AM(t+|At) Q \AtB{t+\At) e \AtB{t+%At) e \AtA(t+\A£) e \AtB{t+\At) _ 

For higher orders we have explicitly: 

%(At) = - 1 -%(At) + ^^y (30) 
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™-^>-£*(t) + £*(t) 

16384^ 4 {At\ 390625 5 f At\ 



2835 z V 4 / 72576 z V 5 / 



7? - +— — Ti - . (33) 



Remark 5 In £/ie worA; of Blanes, Casas and Ros[ll] and Chan and Murua[12] 
the idea of extrapolating symplectic algorithms are also discussed. They pre- 
sented the case of extrapolating an In-order symplectic integrators and noted 
that extrapolating a In-order symplectic integrator will preserve the symplectic 
character of the algorithm to order An + 1 . 



5 Numerical Results 



In the following we deal with the computation of the Hamiltonian, which is 
dervied in Section 2, see also [4,5]. 

Our Hamiltonian of the Levitron is given as: 



H 



-M 



Pl+P2+ P3 



Pi (P5-P6 cosg 4 ) 2 P& 
a a sin 2 g 4 c 



smg 4 coso 5 — smg 5 — + cosg 4 — 
V dqx dq 2 dq 3 



+ <?3 



(34) 



For our splitting scheme, we apply the Hamiltonian of (34), we have: 



OH 
dp 



(p.q) 



(35) 



Pl,P2,P3, 



Pa {Pb -j9ecosg 4 ) 2 Pe(cos 2 g 4 + (a/c) sin 4 ) - p 5 cos g4' 



asm q 4 



asm g 4 
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given as operator A and 



OH 

p = -^(p>q) 



(M sm q 4 cos <?5-?— o- + cos g 4 - — - — , M sin q 4 cos <?5-?— o- + cos g 4 



dqf dqxdqsj \ dq% dq 2 dq 3 

, r , . , <9 2 ^ d 2 ^ \ d 2 #\ 

M sin g 4 smq 5 — — h cos g 5 + cosg 4 — T - 1 

dq 2 dq 3 oqidq-z dq§ ' 



M cosg 4 smg 5 h cosg 5 — - smg 4 — 

V V 92 5 9i/ <9g 3 / 



Pe (P5 - cos g 4 ) cos g 4 (p 5 - p 6 cos g 4 ) 



2 



a sin g 4 a sin 3 q 4 

( ( d^f d^W 
M I sing 4 I cosg 5 — - sing 5 — I I ,0) (36) 

given as operator 5, which we insert into the Algorithm 1. 

We compare our novel schemes (extrapolated Stonier- Verlet method) with 
standard and Runge-Kutta algorithms. Due to the long computation time 
needed, we simulated only 1000 timesteps and compare the trajectory with 
the reference solution from the Runge-Kutta algorithm. In figure 2 is shown 
how the trajectory of the same initial conditions looks like with the Verlet 
algorithm. 




-0.01 -0.008 -0.006 -0.004 -0.002 0.002 0.004 0.006 0.008 



Fig. 2. Trajectory calculated with Verlet algorithm (Left figure: 3D presentation, 
right figure: 2D presentation). 

We improve the solution with an extrapolation scheme with fourth order. 
We have a view at the errors this algorithm produces in comparison with 
the Runge-Kutta Solution with small time-steps (10~ 5 time units per step). 
In Figure 3 and 4, we presented the results of the 4th, 6th and 8th order 
Mult ipro duct expansion method with different time-steps and compared it 
with a fine resolved 4th order Runge-Kutta Benchmark solution (h = 10 -8 ). 

The time scales and computational amount for the extrapolation schemes are 
given in Table 1 and 2. 
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Fig. 3. Errors of the numerical scheme: Extrapolation Scheme with Verlet method 
a Kernel (left figure: 4th order scheme with h = 10~ 5 and 6th order scheme with 
h = IO" 6 ). 





Extrapolation 4th order 


Extrapolation 6th order 


timestep 


1(T 5 


10~ 6 


10~ 5 


10~ 6 


number of steps 


10 s 


10 9 


10 8 


10 9 


computing time 


14min 


142min 


29min 


272min 


mean error 


0.007 


0.007 


0.0068 


0.0068 


maximal error 


0.0226 


0.0234 


0.0188 


0.0188 



Table 1 

Errors and Computational Time with 4th order MPE scheme using Verlet Scheme 
as Kernel. 





Extrapolation 


Extrapolation 


Extrapolation 




6th order 


8th order 


10th order 


timestep 


io- 4 


10~ 3 


io- 2 


10~ 3 


number of steps 


10 7 


10 6 


10 5 


10 6 


computing time 


2.5min 


0.5min 


3sec 


32.8sec 


mean error 


1.0244- 10" 4 


9.6297- 10" 5 


0.013 


2.2397 • IO" 5 


maximal error 


3.4608 • 10" 4 


4.0936 • 10" 4 


0.01 


9.8868 • IO" 5 


le 2 



Errors and Computational Time with higher order MPE scheme using Verlet Scheme 
as Kernel. 

We have also tested the 10th order extrapolation with 10~ 2 time units per 
step and also obtained stable trajectory. 

Remark 6 In the examples, we have verified, that we can improve a basic 
second order symplectic splitting scheme with extrapolation schemes. At least 
achieved higher accurate solutions and save computational time. Moreover we 
save computer resources and obtained stable trajectories with larger time-steps. 
The best result we achieve with the order 10 and h = 10 -2 for such a case we 
could improve the results and are 10-times faster than with standard 4th order 
explicit Runge-Kutta schemes. 



10 



Fig. 4. Errors of the numerical scheme: Extrapolation Scheme with Verlet method 
a Kernel (left figure: 6th order scheme with h = 10 -4 and 8th order scheme with 
h = 1(T 3 ). 

6 Conclusion 



In the paper, we have presented a model to simulate a Levitron. Based on 
the given Hamiltonian system, which can be written as large system of time- 
dependent ordinary differential equation, we present novel and faster solvers 
based on splitting and extrapolation ideas. We could achieve more accurate 
and stable results with higher order schemes and save computational time with 
respect of stable computations. In future, we concentrate on the numerical 
analysis and embedding higher order splitting kernels to nonlinear differential 
equations based on Hamiltonian systems. 
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